/****************************************************************************
*
*					Functions of class SimpleParameterLearning
*
****************************************************************************/

#include <string>
#include <utility>
#include <stdlib.h>
#include <time.h>
#include "GreedyRateRatioSequencingHeuristic.h"

#define eps 0.00000000000000000001

/****************************************************************************
*					GreedyRateRatioSequencingHeuristic::GreedyRateRatioSequencingHeuristic
*
*	Description:	Standard constructor with initialization of internal parameters
*	Parameters:		-
*	Returns:		-
*
****************************************************************************/
GreedyRateRatioSequencingHeuristic::GreedyRateRatioSequencingHeuristic (void) {
	
}


/****************************************************************************
*		GreedyRateRatioSequencingHeuristic::~GreedyRateRatioSequencingHeuristic
*
*	Description:	Standard destructor with removing all dynamic internal variables
*	Parameters:		-
*	Returns:		-
*
****************************************************************************/
GreedyRateRatioSequencingHeuristic::~GreedyRateRatioSequencingHeuristic (void) {
	
}


/****************************************************************************
*					GreedyRateRatioSequencingHeuristic::initStructure
*
*	Description:	Initialize internal variables by observing data set parameters
*	Parameters:		environment - a set of parameters from config file
*	Returns:		0
*
****************************************************************************/
int GreedyRateRatioSequencingHeuristic::initStructure (Environment& environment) {
	double dVal;
	if (environment.getDoubleParamValue(dVal, "gamma") != 0) {
		std::cout << "Error in GreedyRateRatioSequencingHeuristic::initStructure: cannot find 'Gamma' in config file.\n";
		return -1;
	}
	gamma = (float) dVal;
	return 0;
}


/****************************************************************************
*					GreedyRateRatioSequencingHeuristic::findSequence
*
*	Description:	Finds a good (highly probable) sequence, in which the variables 
*					could be changing values between two given timestamps. 
*
*					Also, this function updates the gradients for weights (parameter 
*					Theta) of all ingredients, using formula (for fixed Theta):
*						d(logPseq) / d(Theta) = sum_{events} ( TERM1 - TERM2 / TERM3 ),
*					where 
*						TERM1 = f_Theta(PGMe with i) - f_Theta(PGMe without i),
*					and 
*						TERM2 = sum_{j} lambda_j * (fTheta(PGMe with j) - fTheta(PGMe without j),
*					and
*						TERM3 = sum_{j} lambda_j.
*
*					Since there may exist more than one ingredient, the TERM1, TERM2 and TERM3
*					are recomputed as many times as there are ingredients (or, equivalently, weights).
*					This function works only with variables that can take on two distinct values.
*
*	Parameters:		sequence - a vector of variable ID ordered by the precedence of changes 
*						ensuring the best sequence
*					seqLogLikelihood - the log of the probabiity of the found sequence
*					gradients - a vector of gradients for the ingredients resulting from the found sequence
*					dataSet - an object containing info of at least two data timestamps (additional ones will not be considered) 
*					extendedPGMStruct - an object describing the PGM with variables as determined from and for timestamp0
*	Returns:		0
*
****************************************************************************/
int GreedyRateRatioSequencingHeuristic::findSequence ( 
			std::vector<unsigned int>& sequence, 
			float& seqLogLikelihood, 
			std::vector<float>& gradients,
			DataSet& dataSet, 
			ExtendedPGMStruct& extendedPGMStruct) 
{
	// Debug purposes - calculating time consumed by one run of this procedure
	std::time_t seconds;
	seconds = time (NULL);

	double countEqualCandidates = 1.0;
	unsigned int i, j, bestCandidateInd, bestCandidate = 0;
	std::vector<unsigned int> changes, changes2, vectori;
	unsigned int changesSize = 0, changes2Size = 0;
	vectori.resize(1, 0);
	float rate, score, bestCandidateScore, bestCandidateRate;
	//float TERM1;
	float TERM3;
	bool isX; 
	unsigned int iPos, jPos;
	std::vector <float> ingredientDiffs, TERM2;
	std::map< unsigned int, std::pair<float, std::vector<float>> > revisedRates, bestRevisedRates;
	std::vector < std::pair <unsigned int, float> > revisedRatesVec;
	unsigned int revRatesSize;
	std::map< unsigned int, std::pair<float, std::vector<float>> >::const_iterator it;
	int ret = 0;

	// Build "changes" - a vector of IDs of the variables that *are known* to change from data timestamp 0 to timestamp 1
	// Include the variables that *might* change into vector "changes2"

	std::vector <bool> isInChanges;
	isInChanges.resize(extendedPGMStruct.nVars, false);

	changes.clear();
	changes2.clear();
	unsigned int curVarID;
	unsigned int nActors2 = extendedPGMStruct.nActors * extendedPGMStruct.nActors;
	for (i = 0; i < dataSet.data[0].size(); i++) {
		if (dataSet.data[0][i] != dataSet.data[1][i]) {
			if (i < nActors2) {
				// we found different Y
				extendedPGMStruct.getSquareVarIndex (isX, iPos, jPos, i);
				if (dataSet.data[1][i] == 2) {
					curVarID = extendedPGMStruct.getVarID(false, iPos, jPos);
					changes.push_back(curVarID);
					isInChanges[curVarID] = true;
				}
				else { 
					if (dataSet.data[1][extendedPGMStruct.getSquareVarID(true, iPos, jPos)] == 2) {
						curVarID = extendedPGMStruct.getVarID(false, iPos, jPos);
						changes.push_back(curVarID);
						isInChanges[curVarID] = true;
					}
					else
						changes2.push_back(extendedPGMStruct.getVarID(false, iPos, jPos));
				}
			} else {
				// we found different X
				extendedPGMStruct.getSquareVarIndex (isX, iPos, jPos, i);
				if (iPos < jPos)
					continue;
				if (dataSet.data[1][i] != 0) {
					curVarID = extendedPGMStruct.getVarID(true, iPos, jPos);
					changes.push_back(curVarID);
					isInChanges[curVarID] = true;
				}
			}
		}
	}

	// TEMP PRINTOUT
//	std::cout << "\nVector changes is:";
//	for (unsigned int l = 0; l < changes.size(); l++)
//		std::cout << "\t" << changes[l];
//	std::cout << "\n";

	// If there are no known changes then exit
	if (changes.empty() == true) {
		sequence.clear();
		seqLogLikelihood = 1.0;
		gradients.resize(extendedPGMStruct.weights.size(), 0.0);
		return 0;
	}

	// Revise "changes2" to contain the variables that are not in "changes", and yet, 
	// have the potential to improve the probability of the sequence:
	// The candidates to stay in "changes2" are "Y" variables that have the highest rates in 
	// the timestamp 1

	// Compute rates of variables in "changes2", keep those with rates greater than 1, erase others
	//tempvariables = extendedPGMStruct.variables;
	//extendedPGMStruct.variables = variables2;
	extendedPGMStruct.setVarValuesFromData (dataSet, 1);
	extendedPGMStruct.computeRates();
	i = 0;
	while (i < changes2.size()) {
		if (extendedPGMStruct.rates[changes2[i]].first <= 1.0 ) {
			changes2.erase(changes2.begin() + i);
			continue;
		}
		i++;
	}
	// Make the max number of top elements to remain in "changes2" be no greater than the size of "changes" tripled
	// Remove the "changes2" elements that don't make it into the top
	// Alternatively, let each of the elements in "changes" choose the top 3 candidates that can help them most
	// ( THIS IS NOT DONE YET )
	
	// Change back the values of the "variables", compute the rates for timestamp0
	//extendedPGMStruct.variables = tempvariables;
	extendedPGMStruct.setVarValuesFromData (dataSet, 0);
	extendedPGMStruct.computeRates();


	// Calculating Term2 and Term3
	TERM2.resize(extendedPGMStruct.weights.size(), 0.0);
	for (i = 0; i < extendedPGMStruct.rates.size(); i++) {
		for (j = 0; j < TERM2.size(); j++) {
			TERM2[j] += extendedPGMStruct.rates[i].first * extendedPGMStruct.rates[i].second[j];
		}
	}

	TERM3 = 0.0;
	// Calculate TERM3 (the sum of all the rates after the best candidate has entered the sequence)
	for (i = 0; i < extendedPGMStruct.rates.size(); i++) {
		TERM3 += extendedPGMStruct.rates[i].first;
	}

	// Empty the sequence
	sequence.clear();
	// Clear the gradients
	gradients.resize(extendedPGMStruct.weights.size(), 0.0);

	bestCandidateRate = 0.0; 
	seqLogLikelihood = 0.0;

	// Cycle: for each place in the sequence, until all the required variables are activated
	while (changes.size() > 0) {
		if (changes.size() == 1)
			bool bb =true;
		std::cout << "\tNum of changes left: " << changes.size() << ", ";

		bestCandidate = (unsigned int) -1;
		bestCandidateScore = 0;
		bestCandidateInd = (unsigned int) -1;

		// Cycle: for each variable (changes[i]) in "changes" compute 
		// the score indicating how it increases the overall probability of the elements in "changes"
		for (i = 0; i < changes.size(); i++) {
			// consider a candidate for entering the sequence
			//candidate = changes[i];

			// make sure this variable is not a diagonal link in the adjacency matrix of the data
			extendedPGMStruct.getVarIndex (isX, iPos, jPos, changes[i]);
			if (jPos == iPos) {
				isInChanges[changes[i]] = false;
				changes.erase(changes.begin() + i); // saving checks of the same variables in future
				continue;
			}
			// read the rate of the candidate from the PGM container extendedPGMStruct
			rate = extendedPGMStruct.rates[ changes[i] ].first;
			// if the rate is zero, skip the candidate
			if (rate < eps) {
				//if (changes[i] < 1089)
				//	std::cout << "!";
				continue;
			}

			// begin to assemble the formula for the score of the candidate
			score = rate / (TERM3 - rate);

			// find the revised rates that would result from the candidate entering the sequence
			// the revised rates together with the different in ingredient values
			// are written into map "revisedRates" with candidate as key 
			// and a pair (rate, vector of ingredient differences) as value
			
			//vectori[0] = changes[i];
			//extendedPGMStruct.reviseRates(revisedRates, vectori);
			extendedPGMStruct.reviseRates(revisedRatesVec, revRatesSize, changes[i]);
			//if (revRatesSize > 0) {
				// finish assembling the score of the candidate so that it eventually equals:
				// Score(xi) = Rate(xi) * Prod_{all vars x in vector "change"} (New_Rate(x) + 1) / (Rate(x) + 1) 
				//  * Prod_{all vars y that must not change} (Rate(y) + 1) / (New_Rate(y) + 1)
				//for (it = revisedRates.begin(); it != revisedRates.end(); it++) {
				for (j = 0; j < revRatesSize; j++) {
					// initialize a flag that says whether a variable ("j") is in "changes" or not
					/*isPresent = false;
					for (j = 0; j < changes.size(); j++) {
						if (it->first == changes[j]) {
							isPresent = true;
							break;
						}
					}
					// depending on "isPresent" update the respective product contributing to the score
					if (isPresent) {
						score *= (it->second.first + 1) / (extendedPGMStruct.rates[it->first].first + 1);
					}
					else {
						score *= (extendedPGMStruct.rates[it->first].first + 1) / (it->second.first + 1);
					}
					*/
					/*if (isInChanges[it->first]) {
						score *= (it->second.first + 1) / (extendedPGMStruct.rates[it->first].first + 1);
					} else {
						score *= (extendedPGMStruct.rates[it->first].first + 1) / (it->second.first + 1);
					}*/
					if (isInChanges[revisedRatesVec[j].first]) {
						score *= (revisedRatesVec[j].second + 1) / (extendedPGMStruct.rates[revisedRatesVec[j].first].first + 1);
					} else {
						score *= (extendedPGMStruct.rates[revisedRatesVec[j].first].first + 1) / (revisedRatesVec[j].second + 1);
					}
				}
			//}

			// if the candidate's score beats the score of the previously found best candidate, replace the latter with the former
			if (score > bestCandidateScore + eps) {
				//std::cout << "\t\t\tBetter candidate selected: " << changes[i] << ": " << score << "\n";
				// update and store the best candidate's score, rate, ID, and the changes it would bring to the PGM's "rates" vector
				bestCandidateScore = score;
				bestCandidateRate = rate;
				bestCandidate = changes[i];
				bestCandidateInd = i;
				//bestRevisedRates = revisedRates;
				countEqualCandidates = 1.0;
			} else if (score > bestCandidateScore - eps) {
				// if the candidate's score is equal to the score of the previously found best candidate
				// make the replacement randomly, ensuring that eventually each equal candidate has an equal chance of entering the sequence
				if ( (double) rand() <= RAND_MAX * (1.0 / (1.0 + countEqualCandidates)) ) {
					//std::cout << "\t\t\tEqual candidate selected: " << changes[i] << ": " << score << "\n";
					bestCandidateRate = rate;
					bestCandidate = changes[i];
					bestCandidateInd = i;
					//bestRevisedRates = revisedRates;
					countEqualCandidates++;
				}
			}
		}
		
		// Cycle: for each variable (changes2[i]) in "changes2" see 
		// if it increases the overall probability of the elements in "changes"
		// more than it decreases the probability of the sequence by its own entry 
		// (this cycle mirrors the one above for vector "changes"; the only difference is in the computation of the score)
	/*	for (i = 0; i < changes2.size(); i++) {
			candidate = changes2[i];
			//std::cout << "\t\tCurrent rate: " << i << "\n";
			extendedPGMStruct.getVarIndex (isX, iPos, jPos, candidate);
			if (jPos == iPos)
				continue;
			rate = extendedPGMStruct.rates[candidate].first;
			if (rate < eps)
				continue;

			// There could be used several variants:
			//score = rate / (TERM3 - rate);
			score = rate * rate / ((TERM3 - rate)*(TERM3 - rate));

			// find the revised rates that would result from the change in variable i
			vectori[0] = candidate;
			extendedPGMStruct.reviseRates(revisedRates, vectori);
			if (revisedRates.empty() == false) {
				// compute the score of variable xi:
				// Score(xi) = Rate(xi)^2 * Prod_{all vars x that have to change} (New_Rate(x) / Rate(x)) 
				//  * Prod_{all vars y that must not change} (Rate(y) / New_Rate(y))
				for (it = revisedRates.begin(); it != revisedRates.end(); it++) {
					isPresent = false;
					for (j = 0; j < changes.size(); j++) {
						if (it->first == changes[j]) {
							isPresent = true;
							break;
						}
					}
					if (isPresent) {
						score *= (it->second.first + 1) / (extendedPGMStruct.rates[it->first].first + 1);
					}
					else {
						score *= (extendedPGMStruct.rates[it->first].first + 1) / (it->second.first + 1);
					}
				}
			}

			if (score > bestCandidateScore + eps) {
				//std::cout << "\t\t\tBetter candidate selected: " << candidate << ": " << score << "\n";
				bestCandidateRate = rate;
				bestCandidate = candidate;
				bestCandidateInd = i;
				//bestRevisedRates = revisedRates;
				countEqualCandidates = 1.0;
			}
			else if (score > bestCandidateScore - eps) {
				if ( (double) rand() <= RAND_MAX * (1.0 / (1.0 + countEqualCandidates)) ) {
					//std::cout << "\t\t\tEqual candidate selected: " << candidate << ": " << score << "\n";
					bestCandidateRate = rate;
					bestCandidate = candidate;
					bestCandidateInd = i;
					//bestRevisedRates = revisedRates;
					countEqualCandidates++;
				}
			}
		}
	*/
		// End Cycle

		// printout infromation about the best candidate
		extendedPGMStruct.getVarIndex (isX, iPos, jPos, bestCandidate);
		std::cout << "New element in sequence: ";
		if (isX) std::cout << "X[";
		else std::cout << "Y[";
		std::cout << iPos << "," << jPos << "].\n";

		// update the overall probability of the sequence
		seqLogLikelihood += std::log (bestCandidateRate/TERM3);
		
		// Update gradients (TERM1 and TERM2 are required to be computed for each ingredient)
		for (j = 0; j < gradients.size(); j++) {
			// UPDATE POSSIBLE - update term2 for each j while doing updateRate for bestCandidate
			/*TERM2 = 0.0;
			for (i = 0; i < extendedPGMStruct.rates.size(); i++) {
				TERM2 += extendedPGMStruct.rates[i].first * extendedPGMStruct.rates[i].second[j];
			}*/
			// Calculate TERM1 and TERM2
			//TERM1 = extendedPGMStruct.rates[bestCandidate].second[j];
			if (TERM3 != 0) 
				gradients[j] = gradients[j] + extendedPGMStruct.rates[bestCandidate].second[j] - (TERM2[j] / TERM3);
			else 
				gradients[j] = gradients[j] + extendedPGMStruct.rates[bestCandidate].second[j];
		}

		// Add the variable with the greatest score to the sequence, update the original pgmStruct
		sequence.push_back(bestCandidate);
		TERM3 += extendedPGMStruct.updateRates(bestCandidate, TERM2);

		//extendedPGMStruct.changeVarValue(bestCandidate);
		// Revise the rates in the extendedPGMStruct
	/*	extendedPGMStruct.rates[bestCandidate] = std::make_pair<float, std::vector<float>>(extendedPGMStruct.computeRate(bestCandidate, ingredientDiffs), ingredientDiffs); 
		if (bestRevisedRates.empty() == false) {
			for (it = bestRevisedRates.begin(); it != bestRevisedRates.end(); it++) 
				extendedPGMStruct.rates[it->first] = std::make_pair<float, std::vector<float>>(it->second.first, it->second.second);
		}
	*/
		//extendedPGMStruct.computeRates();

		/* TERM3 = 0.0;
		// Calculate TERM3 (the sum of all the rates after the best candidate has entered the sequence)
		for (i = 0; i < extendedPGMStruct.rates.size(); i++) {
			TERM3 += extendedPGMStruct.rates[i].first;
		}*/

		// If the added variable is in vector "changes", delete it from there 
		if (isInChanges[bestCandidate]) {
			isInChanges[bestCandidate] = false;
			changes.erase(changes.begin() + bestCandidateInd);
		} else {
			// if the best candidate was found in changes2
			changes2.erase(changes2.begin() + bestCandidateInd);
		}
	// End Cycle
	}
	std::cout << "Sequence seqLogLikelihood: " << seqLogLikelihood << "\n";

	// Debug printout
	//std::cout<< "\n\tTime per one run of findSequence is: " << time(NULL) - seconds << " sec.\n";
	
	return 1;
	
}

/****************************************************************************
*					GreedyRateRatioSequencingHeuristic::printoutProgress
*
*	Description:	Function printout current state of sequence identification
*					Specifically it prints current sequence length 
*					and the product of the rates of variables assigned so far
*	Parameters:		-
*	Returns:		-
*
****************************************************************************/
int GreedyRateRatioSequencingHeuristic::printoutProgress(
	void *instance,
    const lbfgsfloatval_t *x,
    const lbfgsfloatval_t *g,
    const lbfgsfloatval_t fx,
    const lbfgsfloatval_t xnorm,
    const lbfgsfloatval_t gnorm,
    const lbfgsfloatval_t step,
    int n,
    int k,
    int ls )
{
#ifdef _DEBUG
	// Report the result
    printf("Iteration %d:\n", k);
    printf("  fx = %f\n", fx);
	for (int i = 0; i < n; i++)
		printf("  x[%d] = %f", i, x[i]);
	printf("\n");
	printf("  xnorm = %f, gnorm = %f, step = %f\n", xnorm, gnorm, step);
#endif
    return 0;
}
